Ước tính Rt của Sởi năm 2024

Serial interval

  • Dữ liệu hiện có không thể tính được Serial interval theo định nghĩa. Do đó, sẽ dựa vào y văn và chọn Serial interval có trung bình là 14.5 ngày và độ lệch chuẩn là 3.25 ngày (Worden et al. 2020).

Ước tính Rt

Code
library(readxl)
library(EpiEstim)
library(ggplot2)
library(dplyr)
library(janitor)
library(tidyr)
library(knitr)
library(lubridate)

df1 <- read_xlsx("data/data.xlsx", sheet = "data")
df1 <- as.data.frame(df1)

Clean data

Code
df1 <- df1 %>% clean_names()

df <- df1 %>% rename(dates = ngay_khoi_phat,
                    ngaysinh = ngay_sinh,
                    ngaynv = ngay_nhap_vien_kham,
                    tiemchung = tinh_trang_tiem_chung,
                    gioi = gioi_tinh)

sum(is.na(df$dates))
[1] 173
Code
# Lấy ngày nhập viện thay cho ngày khởi phát của các ca missing
df$dates[is.na(df$dates)] <- df$ngaynv[is.na(df$dates)]
sum(is.na(df$dates))
[1] 0
  • Có 173 ca bị missing Ngày khởi phát nên sử dụng Ngày nhập viện thay thế cho những ca này.
Code
df$cd <- ifelse(df$phan_loai_chan_doan == "Loại trừ sởi", NA, df$phan_loai_chan_doan)
df <- df[,c("dates", "ngaysinh", "tiemchung", "cd")]

# Loại những ca Sởi loại trừ
df <- na.omit(df)

df_convert <- df %>% group_by(dates) %>% 
  summarise(I = n())

df_complete <- df_convert %>%
    complete(dates = seq(min(dates), max(dates), by = "day")) %>%
    replace_na(list(I = 0))

df_complete$dates <- as.Date(df_complete$dates)

Ước tính hệ số lây nhiễm Rt

  • Bộ dữ liệu bao gồm những ca Sởi xác định (Sởi lâm sàng + Sởi xét nghiệm dương tính) từ ngày 04/03/2024 đến 14/09/2024. Sử dụng ngày khởi phát để ước tính Rt .
Code
mod_raw <- estimate_R(
  incid = df_complete, # Dữ liệu từ 04/03/2024 - end
  method = "parametric_si", 
  config = make_config(
    list(
      mean_si = 14.5, 
      std_si = 3.25
    )
  )
)
Default config will estimate R on weekly sliding windows.
    To change this change the t_start and t_end arguments. 
Warning in estimate_R_func(incid = incid, method = method, si_sample = si_sample, : You're estimating R too early in the epidemic to get the desired
            posterior CV.
Code
df_rt_raw <- mod_raw$R
df_rt_raw$dates <- mod_raw$dates[df_rt_raw$t_end]
df_rt_raw$q1_rt <- df_rt_raw$`Quantile.0.025(R)`
df_rt_raw$q3_rt <- df_rt_raw$`Quantile.0.975(R)`
df_rt_raw$rt <- df_rt_raw$`Mean(R)`
Code
library(patchwork)

p_hist_raw <- ggplot(df_complete, aes(x = dates, y = I)) +
    geom_histogram(stat = "identity", binwidth = 1, width = 1, fill = "lightgrey", color = "black") +
    labs(x = "Day", y = "Incidence") +
    theme_minimal() + 
    theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) + 
    scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-03-04"), ymd("2024-09-14")))

p_raw <- ggplot(df_rt_raw, aes(x = dates)) +
  geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) + 
  geom_line(aes(y = rt), color = "darkorange") +
  geom_hline(yintercept = 1, linetype = "dashed", color = "black") + 
  labs(x = "Day", y = "Rt", title = "04/03/2024 - 14/09/2024") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
  scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-03-04"), ymd("2024-09-14"))) +
  scale_y_continuous(limits = c(0, 20))

p_hist_raw/p_raw

  • Khi ước tính Rt từ 04/03/2024 thì có thể là quá sớm vì các ca bệnh xuất hiện lẻ tẻ từ 04/03/2024 đến khoảng 23/05/2024. Số ca bệnh bắt đầu xuất hiện liên tục từ 23/05/2024 nên sẽ ước tính Rt từ ngày này. Do đó, Rt sẽ được ước tính từ 23/05/2024 đến 10/09/2024.
Code
# Lọc dữ liệu từ 23/05/2024 - end
df_filter <- filter(df_complete, dates >= "2024-05-30")
Code
mod_1w <- estimate_R(
  incid = df_filter, 
  method = "parametric_si", 
  config = make_config(
    list(
      mean_si = 14.5, 
      std_si = 3.25
    )
  )
)
Default config will estimate R on weekly sliding windows.
    To change this change the t_start and t_end arguments. 
Code
# Lấy dữ liệu ra để vẽ Rt
df_rt_1w <- mod_1w$R
df_rt_1w$dates <- mod_1w$dates[df_rt_1w$t_end]
df_rt_1w$q1_rt <- df_rt_1w$`Quantile.0.025(R)`
df_rt_1w$q3_rt <- df_rt_1w$`Quantile.0.975(R)`
df_rt_1w$rt <- df_rt_1w$`Mean(R)`
Code
library(patchwork)

# Biểu đồ đường cong dịch
p_hist <- ggplot(df_filter, aes(x = dates, y = I)) +
    geom_histogram(stat = "identity", binwidth = 1, width = 1, fill = "lightgrey", color = "black") +
    labs(x = "Day", y = "Incidence") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) + 
    scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14")))

# Biểu đồ cho 1 tuần
p_1w <- ggplot(df_rt_1w, aes(x = dates)) +
  geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) + 
  geom_line(aes(y = rt), color = "darkorange") +
  geom_hline(yintercept = 1, linetype = "dashed", color = "black") + 
  labs(x = "Day", y = "Rt", title = "30/05/2024 - 14/09/2024") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
  scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14"))) +
  scale_y_continuous(limits = c(0, 20))
Code
library(patchwork)

(p_hist_raw + p_hist) / (p_raw + p_1w)

  • Trong lệnh esimate_R, chúng ta có thể tùy chọn Sliding window (SW) khác nhau để thể hiện kết quả ước tính Rt. Sliding window là khoảng thời gian chúng ta dùng để xác định Rt, ví dụ nếu chúng ta chọn SW = 7 thì giá trị Rt sẽ được tính bằng trung bình của 7 ngày. Vậy chúng ta nên chọn *SW bằng bao nhiêu? Chúng ta vẽ thử 4 trường hợp của SW lần lượt là: mỗi ngày, mỗi tuần, mỗi 2 tuần và mỗi tháng để lựa chọn.
Code
# SW = 1
t_start <- seq(2, nrow(df_filter)-1)
t_end <- t_start + 1

mod_day <- estimate_R(
  incid = df_filter, 
  method = "parametric_si", 
  config = make_config(
    list(
      mean_si = 14.5, 
      std_si = 3.25,
      t_start = t_start,
      t_end = t_end
    )
  )
)

# Lấy dữ liệu ra để vẽ Rt
df_rt_day <- mod_day$R
df_rt_day$dates <- mod_day$dates[df_rt_day$t_end]
df_rt_day$q1_rt <- df_rt_day$`Quantile.0.025(R)`
df_rt_day$q3_rt <- df_rt_day$`Quantile.0.975(R)`
df_rt_day$rt <- df_rt_day$`Mean(R)`
Code
# SW = 14
t_start <- seq(2, nrow(df_filter)-13)
t_end <- t_start + 13

mod_2w <- estimate_R(
  incid = df_filter, 
  method = "parametric_si", 
  config = make_config(
    list(
      mean_si = 14.5, 
      std_si = 3.25,
      t_start = t_start,
      t_end = t_end
    )
  )
)

# Lấy dữ liệu ra để vẽ Rt
df_rt_2w <- mod_2w$R
df_rt_2w$dates <- mod_2w$dates[df_rt_2w$t_end]
df_rt_2w$q1_rt <- df_rt_2w$`Quantile.0.025(R)`
df_rt_2w$q3_rt <- df_rt_2w$`Quantile.0.975(R)`
df_rt_2w$rt <- df_rt_2w$`Mean(R)`
Code
# SW = 28
t_start <- seq(2, nrow(df_filter)-27)
t_end <- t_start + 27

mod_month <- estimate_R(
  incid = df_filter, 
  method = "parametric_si", 
  config = make_config(
    list(
      mean_si = 14.5, 
      std_si = 3.25,
      t_start = t_start,
      t_end = t_end
    )
  )
)

# Lấy dữ liệu ra để vẽ Rt
df_rt_month <- mod_month$R
df_rt_month$dates <- mod_month$dates[df_rt_month$t_end]
df_rt_month$q1_rt <- df_rt_month$`Quantile.0.025(R)`
df_rt_month$q3_rt <- df_rt_month$`Quantile.0.975(R)`
df_rt_month$rt <- df_rt_month$`Mean(R)`
Code
library(patchwork)

# Biểu đồ đường cong dịch
p_hist <- ggplot(df_filter, aes(x = dates, y = I)) +
    geom_histogram(stat = "identity", binwidth = 1, width = 1, fill = "lightgrey", color = "black") +
    labs(x = "Day", y = "Incidence") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) + 
    scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14")))

# Biểu đồ cho 1 tuần
p_1w <- ggplot(df_rt_1w, aes(x = dates)) +
  geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) + 
  geom_line(aes(y = rt), color = "darkorange") +
  geom_hline(yintercept = 1, linetype = "dashed", color = "black") + 
  labs(x = "Day", y = "Rt", title = "1 Week") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) + 
  scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14"))) +
  scale_y_continuous(limits = c(0, 8))

# Biểu đồ cho 1 ngày
p_day <- ggplot(df_rt_day, aes(x = dates)) +
  geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) + 
  geom_line(aes(y = rt), color = "darkorange") +
  geom_hline(yintercept = 1, linetype = "dashed", color = "black") + 
  labs(x = "Day", y = "Rt", title = "1 Day") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) + 
  scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14"))) +
  scale_y_continuous(limits = c(0, 10))

# Biểu đồ cho 2 tuần
p_2w <- ggplot(df_rt_2w, aes(x = dates)) +
  geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) + 
  geom_line(aes(y = rt), color = "darkorange") +
  geom_hline(yintercept = 1, linetype = "dashed", color = "black") + 
  labs(x = "Day", y = "Rt", title = "2 Week") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) + 
  scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14"))) +
  scale_y_continuous(limits = c(0, 10))

# Biểu đồ cho 1 tháng
p_month <- ggplot(df_rt_month, aes(x = dates)) +
  geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) + 
  geom_line(aes(y = rt), color = "darkorange") +
  geom_hline(yintercept = 1, linetype = "dashed", color = "black") + 
  labs(x = "Day", y = "Rt", title = "1 Month") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
  scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14"))) +
  scale_y_continuous(limits = c(0, 5))

(p_hist + p_hist)/(p_day + p_1w)

Code
library(patchwork)

(p_hist + p_hist)/(p_2w + p_month)

  • Rt cao khi số ca bắt đầu tăng đột ngột và liên tục, càng về sau Rt giảm xuống gần bằng 1 nhưng không duy trì được giá trị Rt < 1 chứng tỏ dịch vẫn đang diễn ra. Với SW = 1 tuần thì Rt ~ 14; Sliding window = 2 tuần thì Rt ~ 19; các khoảng Sliding window khác có giá trị Rt thấp hơn rất nhiều vào ngày ước tính đầu tiên. Theo lý thuyết, Rt sẽ cao khoảng thời gian đầu do dịch mới bùng phát mà 2 tuần trước đó không có ca bệnh xuất hiện liên tục. Do đó, chúng ta chọn giữa SW = 7 hoặc SW = 14 để ước tính.

  • Tuy nhiên, kết quả cho thấy Rt được ước tính với SW là 2 tuần cho kết quả ổn định hơn, giá trị không lên xuống thất thường nên phù hợp để xem xét và đánh giá diễn tích của dịch. Ngoài ra, khoảng tin cậy của Rt trong thời đầu là quá lớn (khi SW là 1 tuần) nên sẽ không đảm bảo tính tin cậy của Rt.

Ước tính Rt với sliding window là 2 tuần

Code
library(patchwork)

p_hist <- ggplot(df_filter, aes(x = dates, y = I)) +
    geom_histogram(stat = "identity", binwidth = 1, width = 1, fill = "lightgrey", color = "black") +
    labs(x = "Day", y = "Incidence") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) + 
    scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14")))

p_2w <- ggplot(df_rt_2w, aes(x = dates)) +
  geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) + 
  geom_line(aes(y = rt), color = "darkorange") +
  geom_hline(yintercept = 1, linetype = "dashed", color = "black") + 
  labs(x = "Day", y = "Rt", title = "13/06 - 14/09") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8)) + 
  scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14"))) +
  scale_y_continuous(breaks = seq(0,10,by = 2))

p_hist/p_2w

  • Dữ liệu để ước tính Rt là từ 30/05/2024, nên Rt sẽ được tính từ ngày 13/06/2024 (Do Serial interval có trung bình là 14 ngày). Rt ~ 7.07 vào 08/06/2024, sau đó giảm dần khi dịch bắt đầu xảy ra (số ca mới xuất hiện liên tục) và cho đến ngày 14/09/2024 (khoảng tuần 36) thì Rt vẫn lớn hơn 1 chứng tỏ dịch vẫn đang diễn ra và số ca bệnh có thể vẫn tiếp tục tăng.
Code
library(knitr)
kable(df_rt_2w[,c("dates", "rt")])
dates rt
2024-06-13 7.0698366
2024-06-14 5.6690111
2024-06-15 4.3508564
2024-06-16 3.2821600
2024-06-17 2.4718480
2024-06-18 2.1445391
2024-06-19 1.8320200
2024-06-20 1.7193623
2024-06-21 1.6256121
2024-06-22 1.5422473
2024-06-23 1.5176277
2024-06-24 1.4954698
2024-06-25 1.3944655
2024-06-26 1.4405630
2024-06-27 1.3418515
2024-06-28 1.1143724
2024-06-29 1.0995479
2024-06-30 1.2085188
2024-07-01 1.4736674
2024-07-02 1.5490903
2024-07-03 1.5791076
2024-07-04 1.5973473
2024-07-05 1.5671345
2024-07-06 1.4948245
2024-07-07 1.4237618
2024-07-08 1.3243248
2024-07-09 1.1656958
2024-07-10 1.0137158
2024-07-11 1.0568080
2024-07-12 1.0603317
2024-07-13 1.1164884
2024-07-14 1.0754648
2024-07-15 0.8953696
2024-07-16 0.9462627
2024-07-17 1.0530632
2024-07-18 1.0882127
2024-07-19 1.1613209
2024-07-20 1.2219629
2024-07-21 1.2947699
2024-07-22 1.4031893
2024-07-23 1.4613128
2024-07-24 1.5730353
2024-07-25 1.6793297
2024-07-26 1.7178643
2024-07-27 1.6554286
2024-07-28 1.6616364
2024-07-29 1.7875178
2024-07-30 1.7556702
2024-07-31 1.6092751
2024-08-01 1.4647907
2024-08-02 1.3967257
2024-08-03 1.3536969
2024-08-04 1.3559618
2024-08-05 1.3996821
2024-08-06 1.5009996
2024-08-07 1.6549297
2024-08-08 1.7120074
2024-08-09 1.8062972
2024-08-10 1.8852197
2024-08-11 1.9690807
2024-08-12 1.8912403
2024-08-13 1.9860504
2024-08-14 1.9307622
2024-08-15 1.9485278
2024-08-16 1.9479556
2024-08-17 1.9839644
2024-08-18 1.9756048
2024-08-19 1.9431634
2024-08-20 1.7682855
2024-08-21 1.6166982
2024-08-22 1.5031903
2024-08-23 1.4052603
2024-08-24 1.3527849
2024-08-25 1.3302986
2024-08-26 1.3423795
2024-08-27 1.2956510
2024-08-28 1.3565252
2024-08-29 1.3408987
2024-08-30 1.2963237
2024-08-31 1.2655504
2024-09-01 1.2486194
2024-09-02 1.2373010
2024-09-03 1.3128349
2024-09-04 1.3294864
2024-09-05 1.3233001
2024-09-06 1.2926711
2024-09-07 1.2373607
2024-09-08 1.1799842
2024-09-09 1.1565229
2024-09-10 1.1126126
2024-09-11 1.0638421
2024-09-12 1.0497415
2024-09-13 1.0363817
2024-09-14 1.0485225